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Abstract 

A linear flux approach is developed for a finite 
element thermal-structural analysis of steady-state 
thermal and structural problems. The element fluxes are 
assumed to vary linearly In the same form as the 
element unknown variables, and the finite element 
matrices are evaluated in closed form. Since numerical 
integration is avoided, significant computational time 
saving is achieved. Solution accuracy and compu- 
tational speed Improvements are demonstrated by 
solving several two- and three-dimensional thermal- 
structural examples. 
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finite element area 

boundary load vector 

material elastic constants 

element matrices, Eq. (6) 

x and y flux components 

convective heat transfer coeff., Eq. (16); 

beam thickness, Eq. (30) 

internal heat generation 

Jacobian matrix 

thermal conductivity 

beam length, Eq. (30) 

components of unit normal vector 

mass matrix 

element interpolation function 
element matrices, Eq. (11) 
heat flux 
load vector 

distance along boundary 
temperature 

fluid recovery temperature 
surface temperature 
reference temperature for zero stress 
surrounding medium temperature 

temperature Increment 
displacement components 
element constants , Eqs. (21)-(28) 
coordinate directions 
coefficient of thermal expansion 
stress components 
Stefan-Boltzmann constant 
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e strain components, Eq. (15); 

emissivity, Eq. (16) 


Subscripts 
T thermal 

s structural 

Superscripts 
T transpose 

Introduction 

For hypersonic vehicles to become a practical 
reality, efficient techniques are needed to analyze light 
weight airframe and engine structures for repeated and 
prolonged exposure to their severe flight environment. 
To understand the structural response under these 
severe aerodynamic loads, research in the Aerothermal 
Loads Branch, NASA Langley Research Center has 
focused on developing effective computational 
approaches for predicting the aerodynamic flow and the 
thermal and structural response of the structure, 
including their interactions^. The approaches consist of 
using: (1) a general automated unstructured gridding to 
discretize the aerodynamic flow field and the structure, 
(2) finite element methods to solve for the environment, 
loads, and response for all three disciplines (flow, 
thermal and structural response), and (3) adaptive mesh 
refinement techniques with error indicators to minimize 
the number of grid points and Increase the solution 
accuracy. 

A Taylor-Galerkin finite element algorithm, has 
been used recently to predict the aerodynamic flow field 
as well as the thermal-structural response for high 
speed flow over leading edges 2 . The approach utilizes: 
(1) a Taylor series expansion in time to establish 
recurrence relations for time marching, and (2) the 
method of weighted residuals with Galerkin's criterion 
for spatial discretization. The governing equations are 
cast in conservation form. The standard primitive 
variables are replaced with their flux counterparts, which 
are assumed to vary linearly over the elements. This 
formulation allows the finite element matrices to be 
evaluated in closed form, thereby avoiding the more 
expensive numerical integration. Since the Taylor- 
Galerkin algorithm Is a time marching (transient) 
algorithm, the full benefits of the linear flux formulation 
for steady-state'problems (steady-state heat transfer and 
static structural problems) has not been exploited. 
Furthermore, most structural problems may be treated as 
quasi static even when the aerodynamic loads and the 
thermal response are transient. 

The purpose of this paper Is to extend the Taylor- 
Galerkin algorithm to steady-state thermal-structural 
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analysis. The linear flux formulation and solution 
procedure are introduced. The finite element matrices 
(in integral forms) which are different from those 
appearing in the conventional finite element formulation, 
are presented. A method to derive these finite element 
matrices in closed form is developed and presented for 
both two- and three-dimensional elements. 

The capability of the linear flux formulation is 
demonstrated using four examples: (1) a thermal 
analysis of a circular plate with internal heat generation, 
(2) a structural analysis of a beam bending due to 
thermal load, (3) a thermal-structural analysis of a two- 
dimensional aerodynamically heated leading edge 
model, and (4) a structural analysis of a three- 
dimensional leading edge model. Results are 
compared with available analytical solutions and the 
conventional finite element solutions. 

Thermal-Structural Formulation 


Solution Procedure 

For simplicity in presenting the linear flux 
algorithm, both the steady-state energy equation and 
structural equilibrium equations are written in the form of 
a scalar equation as, 

|i. + = H (4) 

9x 9y 

Even though the derivation presented below is for the 
thermal analysis, the procedure can be applied directly 
to the structural analysis. 

Linear Flux Assumptions. The key feature of the linear 
flux formulation in the thermal analysis is to assume the 
distribution of the element heat fluxes E and F in the 
same form as the element temperature distribution T, 
that is 


The derivation of finite element equations using a 
linear flux formulation is presented for steady-state 
thermal and structural analyses. For simplicity, the 
derivation presented herein is for two-dimensional 
problems. Extension to three-dimensional problems is 
straightforward. The governing equations are written in 
conservation form so that the linear flux formulation can 
be used directly. This formulation yields finite element 
matrices which can be evaluated in closed form. A 
method to evaluate these closed form matrices is 
described for both two- and three-dimensional elements. 


T(x,y) 

E(x,y) 

F(x,y) 


t N ( x. y ) ] { T } 
[ N ( x. y ) 1 { E } 
t N ( x, y ) ] { F } 


(5) 


where [N(x,y)] are the element interpolation functions, 
and {T}, (E) and {F} are the vectors of the element nodal 
quantities. The assumption of a linear distribution of 
element fluxes E and F which are interpolated in the 
same form as other dependent variables (e.g. T, as 
shown in Eq. (5)), is widely used in the computational 

fluid dynamics 2 . 


Governing Equations 

Heat Transfer. The steady-state thermal response of 
a structure is governed by the energy equation in 
conservation form, 


£< E T> + i7< F T> 


= H, 


(D 


where the subscript T denotes the thermal analysis, Ey 
and Fy are the heat flux components, and Hy is the heat 
source per unit volume. The heat fluxes Ey and Fy are 
related to temperature gradients by Fourier's law. 

Structural Response. The static structural response 
is governed by the equilibrium equations in 
conservation form, 



where the subscript s denotes the structural analysis. 
The vectors {E s } and {F s }, which contain the stress 
components, are given by 


Finite Element Equations. The finite element equations 
are derived using the method of weighted residuals 3 . 
The governing differential equation, Eq. (4), is multiplied 
by the weighting functions, [N(x,y)], and integrated over 
the element area A. Integration by parts is performed to 
produce element integral terms and the boundary 
surface integral terms for application of different types of 
thermal boundary conditions. Details of the derivation 
follows the conventional finite element approach 
described in Ref. 3. The finite element equations 


obtained are in the form, 

[ D x ] { E } + [ Dy ] { F } + {R} + {B} = 0 
In this equation, the matrices [Dx] and [Dy] are 

(6) 


flf)[N]dA 

(7a) 

A 

ID,]- j 

A 


(7b) 

The element nodal vector, [R], associated with the heat 
source, H, is defined as 


{E s } T = K V 

{F,} T = K °y 1 < 3) 

The stress components ox, cy, and txy are related 
to the displacement gradients ana the temperature by 
the generalized Hooke's law. 


{R} = f { N } H dA (8) 

A 

The vector {B} representing the boundary nodal vector is 
defined as 
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{ B) o J { N } [N] ds { I { E } + m{F)) 

8 

- J { N } IN) ds { q} (9) 
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where / and m are the components of a unit vector 
normal to the element boundary. The finite element 
matrices shown In Eqs. (7)-(9) are evaluated In closed 
form as will be demonstrated In the subsequent section. 

Since the fluxes, E and F, are related to the 
temperature gradients given by Fourier's law, 



F - -k f <1°b> 


nodal stress components In two dimensions are 
obtained using the constitutive relations, 

o ( = c (J Ej + P, ( T - T c ) i,j = 1 ,2,3 (15) 

where ej are the nodal strain components, and To is 
the reference temperature for zero thermal stress. The 
material elastic constants, cjj, and the thermal expansion 
parameter, pi, may be temperature dependent. 

Boundary Conditions, The boundary conditions for 
thermal analysis are applied via the boundary nodal 
vector {B} shown in Eq. (9). The vector (q) appearing in 
this equation may be replaced by different types of 
boundary conditions, 

f o (insulated) 


where k is the material thermal conductivity, the element 
nodal flux vectors, (E) and (F), can be expressed In 
terms of element nodal temperature, (T), as 


(E) = - k [Pxl (T) 

(F) = - k [Py] (T) 


(11 a) 
(11b) 


where the matrices, (P x ] and [Py], are related to the 
element shape and are given in the Appendix. The 
element nodal flux vectors, Eqs. (11a) and (11b), are 
substituted into Eq. (6) to obtain the final finite element 
equations In terms of the unknown element nodal 
temperature, (T), In the form, 

[KJ (T) = (R) + (B) (12) 

where the stiffness (or conduction) matrix, [K], Is given by 


q = 


q 9 

h(T s -T f ) 


4 4 . 

E o(T # - T J 


(specified heating) 
(surface convection) 

(surface radiation) 


(16) 


The boundary conditions for the structural 
analysis, such as the applied surface pressure, can be 
added into the structural equations via the surface 
boundary vector. The procedure is Identical to that for 
the thermal analysis previously described and Is 
therefore omitted. 


Derivation of Closed Form Finite Element 
Matrices and Element Nodal Gradients 


[K] = k [Dx] [Px] + k [Dy] [ Py] (13) 

Equation (12) Is In a form similar to that obtained from 
the conventional finite element approach except that the 
latter stiffness matrix Is defined by 

1*1 -| k * 

JklfllflbA I") 

A 

Evaluation of the conventional stiffness matrix, for 
some element types, such as the two-dimensional 
quadrilateral and three-dimensional hexahedral 
elements, requires the use of numerical integration. 

For nonlinear problems (e.g. due to temperature 
dependent thermal conductivity), Eq. (12) Is solved by 
the Newton-Raphson Iteration technique 3 * 4 . This 
procedure Is Identical to that used in the conventional 
finite element approach. 

The approach presented for the thermal analysis 
is applied directly to derive the finite element equations 
for the structural analysis. The equations are Identical to 
Eqs. (5)-(13), where the temperature vector, (T), Is 
replaced by the displacement vector containing the 
components u and v In x- and y-dlrections, and the 
vectors (E) and (F] represent the element nodal stress 
components. For elastic orthotropic materials, typical 


All finite element Integrals, such as [Dx], [Dy], and 
[R), which are given by Eqs. (7) and (8), can be 
expressed in 'closed forms. This is true for simple 
element types (rods, triangles and tetrahedrons) as well 
as for more popular elements (quadrilaterals and 
hexahedrons). In addition, closed form expressions for 
other finite element integrals (such as the consistent 
mass matrix) and the gradients of element variables 
(such as 9T/3x) which are normally required In other 
finite element formulations can be obtained by the 
procedure described below. 

Quadrilateral Element 

Typical finite element integrals for a general 
quadrilateral element, as shown in Fig. 1, are given 
below using natural coordinates 

i i 

[M] =J J [N] t [N] | J | d^dn 07) 

-1 -1 

( D x ] = J. J [N] |J| dn 08) 

•i -i 

where [M] is the consistent mass matrix and [D x ] Is the 
element matrix previously defined In Eq. (7a). The 
determinant of the Jacobian |J| represents the 
transformation from the element global x-y coordinates 
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to the natural coordinates 5-q (see Fig. 1). The 
transformation permits the element Integration to be 
evaluated over a square. The determinant of the 
Jacobian for the two-dimensional quadrilateral element 
is 


N I 


J11 

J12 


J21 

J22 



Jrf.'W 

i,S< N iV 



an example, Vi is obtained by equating Eqs. (22) and 
(19) with the property at node 1 (£=q=-1 In Fig. 1) to give 

Vi = t( x 2 " x l) Y4 + ( X 1 * x 4) Y2 + ( x 4- x 2)Vll /4 ( 23 > 

The use of the determinant of the Jacobian in the 
form of Eq. (21) instead of the original Eq. (19) permits 
the finite element matrices, (Eq. (17) and (18)), to be 
evaluated in closed forms. The use of the symbolic 
manipulation program MACSYMA5 greatly simplified 
this evaluation. 

The determinant of the Jacobian in the form of Eq. 
(21) is also used to derive closed form expressions for 
the element gradients. For example, the element 
temperature gradient, 9T/9x, Is given by 


where 

n, -n,( 5, n)-j (i + §5,)(i+nii|) i=1 4 ( 20 > 

The algebraic expression for the determinant of 
the Jacobian shown in Eq. (19) is in the form of the 
partial derivatives of the element interpolation functions 

Nj(5,q) and the element nodal locations (xj, yi; I = 1 4). 

The expression for the determinant of the Jacobian Is 
quite lengthy (contains a total of 64 terms if fully 
expanded), and thus results in a tedious task for deriving 
the closed form element matrices (Eq. (17) and (18)). 
Such a task becomes almost impossible for the three- 
dimensional 8-node hexahedral element in which the 
determinant of the Jacobian, if fully expanded, contains 
approximately 200,000 terms. To overcome this 
difficulty, the determinant of the Jacobian is rewritten in 
an alternate simpler form, 

| J | = XN^.q) V, (21) 


(19) 


91 

9x 


9N 

, I 3N 

.3$ . 

■ J 12 L ^ . 


) (T) 


(24) 


where J12 arid J22 are defined in Eq. (19). The 
temperature gradient at node 1 can be determined by 

setting £=q=-1 to yield 


9x 1 node 1 


(V r i)V 4 MT r T 4 )y 2 + (T 4 -T 2 ) yi 

(v x i>y 4 + < x i - x 4 )y2 +( v x 2 )y i 


(25) 


where subscripts denote the element node numbers 
shown in Fig. 1. 

The approach presented here is used to derive 
closed form expressions for the other finite element 
matrices (eg. [Dy], and {R}) and nodal variable gradients. 
These closed form expressions are used in the 
formulation for both thermal and structural analyses. 


Hexahedral Element 

The three-dimensional finite element matrices are 
in the same form as shown in Eqs. (17) and (18), except 
that integration is performed over the element volume. 
For an 8-node hexahedral element as shown In Fig. 2, 
the element interpolation functions are, 


where Ni(^.q) are the functions of the natural 
coordinates £ and q, and are selected so that Eq. (21) 
represents the complete order of polynomials of % and 
rj as appearing in the original Eq. (19). The unknown 
constants, Vj, are functions of the nodal coordinates and 
are to be determined. For a two-dimensional 
quadrilateral element, the functions, Nj(^,q), can be 
represented by the Isoparametric quadrilateral element 
interpolation functions, Nj(^,rj), as given by Eq. (20). 
Equation (21) then becomes 


n, = N, ( S/n.Q = (i +^|) 0 + nq,) 0 +K,) 1=1 8 (26) 

and the determinant of the Jacobian is given by, 


Ml 


i(EN,x^ A(EN iy ,) 


|J| . V, 

= N,V, + N 2 V 2 + N3V3 + N 4 V 4 (22) 

The unknowns, Vi , i=1 4, can then be determined 

easily using the properties of Nj(^,q), that is N|(^,q) 
equal to unity at node i and zero at the other nodes. As 


By following the procedure described in the 
previous section, the above determinant of the Jacobian 
Is rewritten In an- alternate simpler form as in Eq. (21). In 
three dimensions, the hexahedral element interpolation 
functions, Nj&q.S), given In Eq. (26). can not be used to 
represent the function, Ni(^,q,Q, in the same fashion as 
In the two-dimensional Jacobian formulation. This is 
because the hexahedral element interpolation functions, 
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Nite.n.C). do not P rovide a complete polynomial a3 
required by the determinant of the Jacobian shown In 
Eq. (27). Therefore the functions, are 

represented by 27-node Langrange cubic element 
interpolation functions to obtain a simpler form for the 
Jacobian determinant, given as Eq. (28) 

1 J | = Z l N,($,ti.C) V, 1=1 27 (28) 


analytical method, linear flux approach, and the 
conventional approach are shown in Table 1. These 
results indicate that the linear flux formulation provides 
slightly higher solution accuracy compared to the 
conventional finite element formulation. This Is due to 
the fact that four Gauss point numerical Integration, 
which is commonly used to evaluate the conventlona 
finite element stiffness matrix, (see Eq. (14)), can not 
provide exact Integration for arbitrary quadrilateral 

clamant ehnnPS. 


The unknowns Vj are determined by equating the Eqs. 
(27) and (28) at the nodal locations of the Lagrange 
cubic element. The program MACSYMA 6 was used to 
derive these unknowns, Vj , as well as the associated 
element matrices and the element nodal gradients in 
closed form. The algebraic expressions for the matrices 
and the gradients were very lengthy therefore they were 
translated into FORTRAN statements using MACSYMA 
and were used directly in the analysis code. The use of 
these closed form expressions reduces the 
computational time compared to the traditional 
numerical integration, as will be demonstrated in the 
next section. 

Application 

Four examples are presented to demonstrate the 
accuracy and computational efficiency of the linear flux 
formulation. These consist of: (1) a thermal analysis of a 
circular plate with internal heat generation, (2) a 
structural analysis of a beam bending due to thermal 
load, (3) a two-dimensional thermal-structural analysis 
of an aerodynamically heated leading edge model, and 
(4) a three-dimensional structural analysis of an 
aerodynamically heated leading edge model. Results 
obtained by the linear flux formulation are compared 
with available analytical solutions and the conventional 
finite element solutions. 

Circular Plate With Internal Heat Generation 

A 20 in. diameter stainless steel circular plate with 
internal heat generation and specified zero temperature 
along the circumferential boundary shown in Fig. 3 is 
used as the first example. Analytical solution for the 
variation of temperature in the plate is available and is 
given by 

T(r) = -^r (100-(* ) (29) 

where Q is uniformly distributed internal heat generation 
■ rate per unit volume, k the plate thermal conductivity and 
r the plate radial distance. Due to symmetry, only a 
quarter of the plate is modeled with 10 quadrilatera 
elements as shown in Fig. 3. The distorted quadrilateral 
element shape was selected to evaluate the formulation 
performance under an arbitrary unstructured mesh 

conditio^ Q temperature distribution obtained from 
the linear flux formulation Is compared with the 
analytical solution (Eq. (29)) and the conventional finite 
element solution In Fig. 4. Because of the crude mesh 
used in the model, both finite element solutions 
underpredict the temperature distribution. The tempe- 
ratures at nodes A, B, and C (see Fig. 3), obtained by the 


Table 1 Comparative nodal temperature (°F) 
and errors 


Location 

Analytical 
Eq. (29) 

Linear flux 
Temp. Error 
% 

Conventional 
Temp. Error 
% 

A 

100.0 

96.9 

3.1 

96.5 

3.5 

B 

92.0 

88.3 

4.0 

87.9 

4.4 

C 

68.0 

67.8 

0.3 

67.3 

1.0 


The computational time for linear flux approach 
and the conventional approach are given In Table 2. 
The comparison of the computational time Indicates a 
39% savings for the linear flux approach. The time 
savings is due to the use of the closed form algebraic 
expressions rather than the numerical integration to 
evaluate the finite element stiffness matrix. 

Table 2 Comparative CPU time (CRAY-2 
seconds) for evaluating element 
stiffness matrix 


Linear flux 

Conventional 

% Saving 

0.1036 x 10 * 3 

0.1690x10-3 

39 


Beam Bending Due to Thermal Load 


As a second example, a 4 In. long, 0.1 in. thick 
stainless steel beam pinned on the bottom edges s 
considered and is shown In Fig. 5. The beam Is 
assumed to be flat, and stress free at room temperature. 
The beam temperature Is raised uniformly by 65 F. The 
edge constraints cause the beam to bend Into a convex 
shape. At this relatively low temperature and small 
deformation, beam structural response may be 
approximated by the beam-column theory 6 , In which 
shear effects are neglected. For cylindrical bending, 
flexural rigidity of the beam, D, Is equal to Eh 3 /12, where 
E Is the modulus of elasticity and h Is the thickness of the 
beam (see Fig. 5). The deflection, v(x), Is given by 7 

v(x)= Jj-(«>sM*1) (30) 

2 cosf 

where I Is length of the beam and \ Is V (P/D) 1 . The axial 
constraint force, P, is computed from 
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Fih^n'iy 

/ sin XK 


]£ P_ [ 2 ' 2X 
® D cos 2 ^ 


■) 


- a AT l 


(31) 


where AT Is the beam temperature Increment. 

Using symmetry, one-half of the beam Is modeled 
with 160 quadrilateral finite elements. The finite 
elements are uniformly distributed with 40 elements 
along the beam length and 4 elements through the 
beam thickness. Both the linear flux . and he 
conventional finite element approaches yield Identical 
beam deflections. The predicted deflection distribution 
Is compared with the solution from the beam theory, (Eq. 
(30)), In Fig. 6. The figure shows a very good agreement 
of the beam deflection distributions with the maximum 
difference of about 2% at the beam center (x=0). 

The predicted beam deflection obtained from the 
linear flux formulation demonstrates the capability of the 
approach for providing the same solution accuracy as 
the conventional finite element approach (or the 
structural analysis. Of course, the computational time 
saving of 39% Is still achieved in the evaluation of the 
finite element stiffness matrix. 

Two-Dimensional Leading Edge Model 

To further demonstrate the capability of the linear 
flux formulation for both thermal and structural analyses, 
the approach is applied to predict the thermal-structural 
response of an aerodynamically heated leading edge 
subjected to a high speed flow. The leading edge 
consists of a 0.25 In. nose diameter, 3 In. long model 
made of 0.1 in. thick Inconel 617 alloy as shown In Fig. 
7 The thermal boundary conditions along the outer 
surface consist of applied aerodynamic heating and 
emitted surface radiation. The leading edge Is Insulated 
along the inner surface. A schematic of the finite 
element model, which consisted of 508 quadrilateral 
elements is given In the figure. The mesh Is graded with 
five elements through the thickness and 127 elements 
along the circumference. Approximately 70% of the 
elements lie in the 0.25 in. nose of the leading edge. 

The aerodynamic heating along the leading edge 
outer surface and the aerodynamic flow field 
represented by the Mach number contours are shown in 
Fig 8. These flow solutions were obtained from Ref. 8 
by solving the Navler-Stokes equations. The Mach 
number contours Indicate an unsymmetrlc bow shock 
shape from the free stream Mach 5.25 flow, which Is 
Inclined 12.5° relative to the bottom of the leading edge. 
The aerodynamic heating rate distribution is relatively 
low along both lower and upper surfaces of the leading 
edge compared to the stagnation point heating rate. 
The aerodynamic heating rate on the lower surface is 
slightly higher than the upper surface. The aerodynamic 
heating Increases significantly at the nose because the 
flow stagnates In that region. , . 

The leading edge aerodynamic heating, shown in 
the Fig. 8, was predicted assuming a uniform surface 
temperature of 530°R. During the transient response, 
the aerodynamic heating rate decreases as the leading 
edge temperature Increases. Thus, to obtain a rea Istlc 
leading edge temperature response, the specified 


aerodynamic heating (qs in Eq. (16)) is converted Into 
the surface convection boundary condition (h(Ts-Tr) in 
Eq. (16)). It should also be noted that the change In the 
surface convection coefficient, h, with the surface 
temperature, T s , is small compared to the change In the 

heating rate, qs- , 

The predicted steady-state leading edge 
temperature contours and the outer surface temperature 
distribution are compared with the conventional finite 
element solutions in Fig. 9. The temperature 
distributions obtained from both approaches are almost 
identical with the maximum difference of 0.2% at the 
nose of the leading edge where the peak temperature 

occurs. , , ., . 

The aerodynamic pressure on the leading edge 

and the flow pressure contours® are shown in Fig 10. 
The peak pressure occurs at the flow stagnation point on 
the nose and the pressure is nearly uniform on the top 
and bottom surfaces. This aerodynamic pressure (Fig. 
10) and the leading edge temperature (Fig. 9) are used 
as the aerothermal loads for prediction of the leading 
edge structural response. The leading edge material 
properties, such as the modulus of elasticity and the 
coefficient of thermal expansion, are temperature 
dependent 9 . The finite element discretization previously 
used In the thermal analysis is also used for the 
structural analysis to eliminate the data manipulation 
normally required by the different analysis disciplines. 
The predicted tangential stress contours superimposed 
on the deformed leading edge are shown in Fig. 11. 
The increased leading edge temperature causes the 
leading edge to expand. The temperature difference 
between the lower and upper sections (higher 
temperature on the lower section) causes the leading 
edge to bend and rotate upward. The figure shows that 
the tangential stress, which Is primarily caused by the 
temperature difference between the two sections, is 
relatively low. These results are In excellent agreement 
with the conventional finite element results Indicating the 
validity of the’ linear flux formulation for the structural 
analysis. Again, a 39% computational time saving is 
achieved using the linear flux approach. 

The predicted structural response obtained from 
the two-dimensional leading edge model is based on a 
plane strain assumption. The use of this assumption 
results in high compressive axial stresses (-150 ksi) in 
the direction normal to leading edge cross-section. A 
three-dimensional analysis with the appropriate 
boundary conditions would provide a more realistic 
leading edge axial stress prediction. Such a structural 
analysis is presented In the next application. 

Three-Dimensional Leading Edge Model 

As mentioned in the theoretical formulation 
section, the extension of the linear flux approach to three 
dimensions Is straightforward. The approach has been 
extended for both the thermal and structural analyses 
using the 8-node hexahedral element. The use of the 
hexahedral element Is preferred over other three- 
dimensional elfement types (such as the tetrahedral 
element) to reduce the computer memory required for 
the analysis (a hexahedral element consist of five 
tetrahedral elements). The purpose of presenting this 
application is: (1) to compare the linear flux solution with 
the conventional finite element solution for a three- 
dimensional problem, (2) to demonstrate the 
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computational time saving, and (3) to predict a more 
realistic leading edge axial stress. 

The linear flux approach is applied to predict the 
leading edge structural response using the three- 
dimensional model shown in Fig. 12. The mesh on the 
leading edge cross-section (x-y plane) is identical to the 
two-dimensional model described in the previous 
example. The mesh is extended with a total of 10 layers 
in the z-direction. The aerothermal loads consist of the 
temperature distribution (Fig. 9) and the aerodynamic 
pressure (Fig. 10) which are taken to be uniform in the z- 

direction. . 

The predicted axial stress contours superimposed 
on the deformed leading edge along the midsection 
(section z=0.75 in.) are shown in Fig. 13. This figure 
shows a more realistic axial stress distribution which 
resembles the temperature distribution. The peak 
compressive stress of approximately 20 ksi (compared 
to 150 ksi for 2-D model) occurs at the nose of the 
leading edge where the temperature is maximum. This 
significant reduction occurs because there is no axial 
constraint as the two-dimensional case. The 
conventional finite element analysis solution is nearly 
identical with a maximum difference in the deflection at 
the leading edge nose of less than 0.5%. The linear flux 
approach shows a computational time saving of 78% 
over the conventional finite element approach in the 
evaluation of the hexahedral element stiffness matrix. 
Such a significant computational time saving is due to 
the use of closed form algebraic expressions instead of 
the 8-Gauss point numerical Integration. Of course, 
fewer Gauss point integration could be used at the 
expense of accuracy. 

Concluding Remarks 

A linear flux approach for finite element thermal- 
structural analysis was presented. The approach 
employs the assumption that the fluxes vary linearly as 
the dependent variables over the element.. Such an 
assumption is widely used in computational fluid 
dynamics. The finite element equations for steady-state 
thermal and structural analyses are derived. The finite 
element equations consist of the finite element matrices 
in integral form which are different from those appearing 
in the conventional finite element formulation. A method 
was developed to derive these finite element matrices in 
closed form and the details of the derivation was 
described. The use of the closed form algebraic 
expressions for evaluating the finite element matrices 
reduces the computational time by 39% for the 2D 
problems and 78% for the 3D problems compared to 
numerical integration. The linear flux formulation also 
yields slightly higher accurate results compared to the 
conventional finite element formulation. 

Four thermal and structural problems were 
analyzed, and the results compare favorably with 
available analytical and conventional finite element 
results. The examples have demonstrated the viability 
of the approach to improve the disciplinary analysis 
efficiency for practical steady-state thermal-structural 
problems. 


Appendix 

Closed Form Finite Element Matrices 


The closed form expressions of the finite element 
matrices for two-dimensional quadrilateral element, 
such as the [D x ] matrix, shown In Eq. (7), are 


Dx ( 1 . 1 ) 

Dx ( 2 , 2 ) 

Dx ( 1 ,3) 

Dx ( 2 , 4 ) 

Dx(1 .2) 

Dx ( '1 . 4 ) 
Dx ( 2 , 1 ) 
Dx ( 2 , 3 ) 
Dx ( 3 . 2 ) 

Dx ( 3 , 4 ) 
Dx ( 4 , 1 ) 
Dx ( 4 , 3 ) 


= 

-D x (3.3) 

= - 

(y 4 -y 2 )/6 

= 

* D x ( 4 , 4 ) 

= - 

(yi - y 3 ) / 6 

= 

~ D x (3,1) 

= - 

(y4 - y2 ) / 1 2 

= 

-D x (4,2) 

= - 

(yi-y 3 )/12 


- ( y4 + 

y3 

- 2 y 2 ) / 1 2 

= 

- (2 y 4 - 

ya 

- yz)/i2 

= 

( y4 + 

ya 

- 2 y,)/12 

s 

- ( y4 - 

2 y 3 

+ yi)/i2 

= 

( y4 • 

2 y 2 

- yi)/i2 

= 

(2 n • 

ya 

- yi)/i2 

= 

- ( n + 

yz 

- 2 yi)/12 

rs 

- (2 y 3 - 

Vz 

- yi)/i2 


where xj and yj, i=i 4 are nodal coordinates of 

the element based on the element node numbering 
scheme shown In Fig. 1. 

The element flux variation and the element nodal 
variables are related through the matrices [Px] and [Py] 
given in Eq. (11). The matrix, [P x ], as an example, is 
given by 


p,(i.i) = (y 2 -y 4 ) /4V i 

p x (i.2) = (y 4 -y i )/4v 1 
p x (i .4) = (y^yjOMV, 
p x ( 2 , i ) = (y 2 -y 3 )/4V 2 

p x (2,2) = (y 3 -y 1 )/4V 2 
p x ( 2 , 3 ) = (y 1 -y 2 )/4V 2 
p x ( 3 . 2 ) = (y 3 -y 4 )/4V 3 
p x (3,3) = (y 4 - y 2 ) / 4 v 3 
P x ( 3 , 4 ) = (y 2 -y 3 )/4V 3 
p x ( 4 1 1 ) = (y 3 -y 4 )/4V 4 
P x (4,3) = (y 4 -yi)/4V 4 
Px(4.4) = (y r y 3 )/4V 4 
P x (1 ,3) = Px(2,4) = P x (3,1 ) = P X (4,2) = ° 

where V| ( |=i 4 are the constants shown In Eq. (22), 

and are given by, 

Vi = [(x z -x 1 )y 4 + (x,-x 4 )y 2 + (x 4 -x 2 )y 1 ]/4 

V 2 = l( x 2 -x 1 )y 3 + (x, - x 3 ) y 2 + (x 3 - x 2 ) y,]/4 

V 3 = I(x 3 - x 2 ) y 4 + (x 2 - x 4 ) y 3 + (x 4 -x 3 )y 2 ]/4 

V 4 = [(x 3 -x 1 )y 4 + (x, - x 4 ) y 3 + (x 4 -x 3 )y 1 ]/4 


7 



References 


1 . Wieting, A. R. , Dechaumphai, P., Bey, K. S., 
Thornton, E. A. and Morgan, K., "Application of 
Integrated Fluid-Thermal-Structural Analysis 
Method," Presented at the 16th Congress of the 
International Council of the Aeronautical Sciences, 
Jerusalem, Israel, August 28- September 2, 1988. 

2. Dechaumphai, P., Thornton, E. A. and Wieting, A. 
R„ "Flow-Thermal-Structural Study of 
Aerodynamically Heated Leading Edges," AIAA 
Paper No. 88-2245-CP, April 1988. 

3. Huebner, K. H. and Thornton, E. A., The Finite 
Element Method for Engineers, Wiley, New York, 
1982. 

4. Reddy, J. N., An Introduction to the Finite Element 
Method, McGraw-Hill, New York, 1984. 



Fig. 2 Hexahedral finite element In three-dimensional 
global and natural coordinates. 


5. The Mathlab Group, MACSYMA Reference 
Manual, The Mathlab Group Laboratory for 
Computer Science, MIT, Version 9, Second 
Printing, December 1977. 

6. Boley, B. A. and Weiner, J. H., Theory of Thermal 
Stresses, Wiley, New York, 1960. 

7. Thornton, E. A. and Dechaumphai, P., "Coupled 
Flow, Thermal, and Structural Analysis of 
Aerodynamically Heated Panels," Journal of 
Aircraft, Vol. 25, No. 11, November 1988, pp. 1052- 
1059. 

8. Dechaumphai, P., Wieting, A. R. and Pandey, A. K., 
"Flow-Thermal-Structural Interaction of 
Aerodynamically Heated Leading Edges," AIAA 
Paper No. 89-1227-CP, April 1989. 

9. Metals Handbook Committee, Metals Handbook, 
Eighth Edition, American Society for Metals, Ohio, 
1975. 



Fig. 3 Finite element model for a 10-Inch radius 
circular plate subjected to Internal heat 
generation. 



Fig. 4 Comparative temperature distribution for 1 0- 
Inch radius circular plate subjected to Internal 
heat generation. 


Fig. 1 Quadrilateral finite element In two-dimensional 
global and natural coordinates. 


8 




160 elements 
y 210 nodes 



Fig. 5 A schematic finite element model of heated 
beam with boundary conditions. 
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Fig. 8 Surface heating rate distribution and flow Mach 
number contours for an undisturbed Mach 5.25 
flow over leading edge. 
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Fig. 6 Comparative deflections for a 4 in. long beam Fig. 9 Steady-state surface temperature distributions 
subjected to uniform temperature Increment of and leading edge temperature contours. 

65°F. 



s, In. 


Fig. 7 A schematic thermal-structural finite element Fig. 10 Surface pressure distribution and flow pressure 
model of 0.25-Inch diameter, 3-Inch long contours for an undisturbed Mach 5.25 flow over 

leading edge with boundary conditions. leading edge. 
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Fig. 1 1 Tangential stress contours on deformed leading 
edge. 



Fig. 12 A schematic three-dimensional finite element 
model of 0.25-inch diameter, 3-inch long, 
1.5-inch wide leading edge with boundary 
conditions. 



Fig. 13 Axial stress contours on deformed leading 
edge. 
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